#include "solver_stable.h"

#include <math.h>

double* solve_stable(double* u0, Area area, double epsilon, TwoDimFunc f)
{
    int i, j;
    int n = area.n;
    double a = area.a, b = area.b;
    double h = (b - a) / n;
    double dmax, d, dm, temp;

    double* u = new double[n * n];
    double* un = new double[n * n];

    for (i = 0; i < n; ++i)
    {
        for (j = 0; j < n; ++j)
        {
            u[i * n + j] = u0[i * n + j];
        }
    }

    int k = 0;
    do 
    {
        dmax = 0;

        for (i = 1; i < n - 1; i++) 
        {
            dm = 0;
            for (j = 1; j < n - 1; j++) 
            {
                temp = u[i * n + j];
                un[i * n + j] = 0.25 * (
                    u[(i - 1) * n + j] +
                    u[(i + 1) * n + j] +
                    u[i * n + j - 1] +
                    u[i * n + j + 1] -
                    h * h * f(a + i * h, a + j * h));
                d = abs(temp - un[i * n + j]);
                if (dm < d) dm = d;
            }
            if ( dmax < dm ) dmax = dm;
        } 

        for (i = 1; i < n - 1; i++)
        {
            for (j = 1; j < n - 1; j++)
            {
                u[i * n + j] = un[i * n + j];
            }
        }
    } while (dmax > epsilon);

    delete[] un;
    return u;
}